Accepted for publication in ApJ 

Resolving the Circumstellar Disk of HL Tauri at Millimeter 

Wavelengths 

Woojin Kwon and Leslie W. Looney 

^ ! Department of Astronomy, University of Illinois, 1002 West Green Street, Urbana, IL 



(N 

00 



o 



61801 

wkwon@illinois . edu 

and 

Lee G. Mundy 



F^ . Department of Astronomy, University of Maryland, College Park, MD 20742 



ABSTRACT 



We present results of high-resolution imaging toward HL Tau by the Com- 
bined Array for Research in Millimeter- wave Astronomy (CARMA). We have 
. ^ ! obtained A = 1.3 mm and 2.7 mm dust continua with an angular resolution down 

f^>. , to 0.13". Through model fitting to the two wavelength data simultaneously in 

^ ! Bayesian inference using a flared viscous accretion disk model, we estimate the 

physical properties of HL Tau, such as density distribution, dust opacity spectral 
O I index, disk mass, disk size, inclination angle, position angle, and disk thickness. 

HL Tau has a circumstellar disk mass of 0.13 Mq, a characteristic radius of 79 
AU, an inclination of 40°, and a position angle of 136°. Although a thin disk 
k/< I model is preferred by our two wavelength data, a thick disk model is needed to 

H I explain the high mid- and far-infrared emission of the HL Tau spectral energy 

distribution. This could imply large dust grains settled down on the mid plane 
with fine dust grains mixed with gas. The HL Tau disk is likely gravitationally 
unstable and can be fragmented between 50 and 100 AU of radius. However, 
we did not detect dust thermal continuum supporting the protoplanet candi- 
date claimed by a previous study using observations of the Very Large Array at 
A = 1.3 cm. 

Subject headings: circumstellar matter — planetary systems: protoplanetary 
disks — radio continuum: stars — stars: individual (HL Tau) — stars: pre- 
main-sequence — techniques: interferometric 
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INTRODUCTION 



Young circumstellar disks are observed over a wide range of wavelengths to study the 
physical conditions in the disk, the distribution of the gas and dust, and changes in the prop- 
erties of the disk material. The thermal emission from dust at millimeter to sub-millimeter 
wavelengths is the best probe of the bulk material distribution in the disk and the best 
monitor of the growth of large grains within the disk. Due to the high angular resolutions 
available with radio interferometers, continuum observations at these wavelengths provide 
some of the best constraints on the density distribution and the overall st ructure of disks 



durin g the early era when giant planet formation is likely occurring (e.g., iD'Angelo et al. 
2010h . 



This paper presents sub-arcsecond imaging of the disk around HL Tau at millimeter 
wavelengths with the Combined Array for Research in Millimeter- wave Astronomy (CARMA). 
HL Tau is a T Tauri star in the Taurus molecu lar cloud, a nearby star forming region 
at a distance of 140 pc (e.g., iRebuU et al.l |2004[ ) . HL Tau has been studied extensively 



In- 



over the last two decades. The disk has been studied by observations in sub/millimeter 
wavelength continuum using various single dishes and interferome ters: for examp 
stitut de Radioastronomie Millimetrique (IRAM) at A = L3 mm (IBeckwith et al. 
Berkeley-IUinois-Ma ryland Association array (BIMA) at A = 2.7 mm ( jMundy et al. 



199(1) . 



11996 



Loonev et al. 120001) , Very Large Array ( VL A) at centimeter and millimeter wavelengths 
( IWilner et al.l Il996t iGreaves et al.l l2008l ). Owens Valley Radio Observatory (OVRO) and 
interferometry of Caltech Submillimeter Obser yatory (CSO) a nd James Clerk Maxwell Tele- 



scope (JCMT ) at sub /millimeter wavel engths (JLay et a 
wavelengths (jChandler &: Richerll2000l ). In particular. 



1997). and JCMT at submillimeter 



Mundy et al.l ( 11996! ) found physical 



properties such as density and temperature distribution, disk mass, outer radius, and incli- 
nation angle, through power-law disk model fitting to sub-a rcsecond (1'.'32 x 0'.'48) angular 
resolution data of BIMA at A = 2.7 mm, and JLay et al.l (119971) studied t he HL Tau disk prop- 
erties using Bayesian inference for the first time. iKitamura et al.l ( 12002| ) also studied HL Tau 
by modeling of a viscous accretion disk as well as a power-law disk. In addition to the dis k 
structure, HL Tau has an optical jet and molecular bipolar outflow (e.g.. lMundt et al.lll990l ). 



Welch et al.l (l2000l ) showed that HL Tau is located in a wall of a bubb le structure usin g CO 



observations of the BIMA array and the NRAO 12 m telescope and JRobitaille et al.l (120071 ) 
reported that HL Tau has an envelope component through modeling to spectral energy 
distribution over optical, infrared, and submillimeter wavelengths. 



Greaves et al.l ( l2008l ) claimed to detect a protoplanet candidate at a projected radius of 
55 AU using VLA at A = 1.3 cm with 0.08" ang ular resolution; the position is coincident with 
a secondary peak in previous A = 1.3 mm data (jWelch et al.ll2004l ). JNero fc BjorkmanI (120091 ) 
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discussed the possibility of planet formation through disk fragmentation where iGreaves et al. 
( 120081 ) found a compact feature in HL Tau. Based on the Toomre Q parameter and perturba- 
tion cooling time, they argue that the c ompa ct feature could be a result of disk fragmentation. 
In contrast, ICarrasco-Gonzalez et al.l (120091 ) did not detect an emission peak corresponding 
to the protoplanet candidate, at A = 7 mm with ~ 0.05" angular resolution. 

In this paper, we map the structure of the disk around HL Tau and reveal its properties 
by visibility modeling. We briefly discuss how the CARMA data have been taken and reduced 
in Section |21 and our disk modeling is presented in Section |3J Observational and modeling 
results, and the implications for the candidate protoplanet are described and discussed in 
Section HI followed by conclusion in Section [51 



OBSERVATIONS AND DATA 



We carried out observations toward HL Tau in the A = 1.3 mm continuum using the A, 
B, and C configurations and in the 2.7 mm continuum using the B and C configurations of 
CARMA JWoodvet al.ll2004h . The datasets were taken between 2007 November and 2009 
January using the 10.4-m and 6.1-m antennas. To check the gain calibration, particularly for 
extended-array data of A and B configurations, a test calibrator is employed to verify success- 
ful calibration. In add ition, we used the CARMA Paired A ntenna Calibration System for the 
A configuration data (JLamb et al.ll2009l : iPerez et al.ll2010l ): the 3.5-m antennas continuously 
observe a calibrator and their data are used to correct short atmospheric perturbations. The 
improvement of calibration was about 10-20% in terms of image noise levels and the size, 
flux, and peak intensity of the test calibrator. Telescope pointing during th e observations 
were monitored using optical pointing and radio pointing ( ICorder et al.l 120101 ). To minimize 
the bias induced by flux calibration uncertainty, a common flux calibrator (e.g., Uranus) 
was used to bootstrap gain calibrator fluxes. In addition, different array-configuration data 
for HL Tau were compared at common uv distances. Sub/millimeter dust emission from T 
Tauri disks is not variable over a period of a few years, so amplitudes at common uv positions 
should be comparable even in different configuration arrays. We estimate that the absolute 
flux calibration uncertainty is 10% at A = 1.3 mm and 8% at A = 2.7 mm. 



MIRIAD (jSault et al.lll995l ) was employed to calibrate and map data. In addition to 
normal calibration, seeing has been corrected using the UVCAL task for the B configu- 
ration data at A = 1.3 mm, which were taken under less favorable weather conditions. 
Individual configuration data were calibrated separately and combined to make maps. The 
proper motion of HL Tau {vra = 8.0 ± 6.0 mas year~^ and voec = —21.8 ± 5.8 mas year~^; 
Zacharias et al.ll2003l ) was compensated for the data to set the epoch positions to 2009 Jan- 
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uary. The sensitivity and emphasized size scales in the maps depend on weighting schemes of 
visib ihty data. In order to emphasize the small structures, we used Briggs robust weighting 
of teriggslllQQsl ). 



Figure [T] shows the A, B, and C configuration-combined image at A = 1.3 mm (z/ = 229 
GHz) with O'.'ir X O'.'IS (PA = 85.0°) resolution, 24 AU by 18 AU in linear size at the 
distance of HL Tau. The RMS noise level in the A = 1.3 mm map is 0.8 mJy beam~^. The 
HL Tau disk is nicely resolved; the apparent disk size is ~ I'.'S x I'.'l at the Aa contour 
level, corresponding to 210 AU major axis. The peak intensity at A = 1.3 mm is 33 mJy 
beam~^ and the total flux in a 2" box centered on the source is 700 ± 10.3 mJy. A Gaussian 
fit to the emission yields an emission centroid position of RA(J2000) = 04^^ 31^" 38!418, 
Dec(J2000) = +18° 13' 5r.'37. The position angle of the major axis is 135° East-of-North 
and the inclination is 43°, based on the ratio of the major and minor axes (0° corresponds 
to a face-on disk). The A = 2.7 mm (z/ = 112 GHz) data yield a continuum map with a 
resolution of I'.'Ol x 0'.'67 (PA = 72.8°) and an RMS noise level of 1.1 mJy beam~^; the peak 
intensity is 64 mJy beam~^ and the integrated flux is 120 ± 3.8 mJy. The given errors in 
both fluxes are statistical uncertainties; systematic calibration uncertainties are estimated 
to be about 10%. 



3. DISK MODELING 



The disk around HL Tau shows a high degree of axial symmetry and a significant extent 
relative to the beam size at A = 1.3 mm wavelength. In order to derive physical para meters 
for the disk, we fit the data with a standard viscous accretion disk model (jPringldl 19811 ). This 
disk model has a power-law radial density distribution tapered by an exponential function: 
in the case of a thin disk, S(-R) oc {R/ Rc)~'^exp[—{R/ RcY~^], where Rr. is a characteristic 



radius (e.g., [Andrews et al.ll2009l ). For our fits, we assume a disk thickness determined by 



vertical hydrostatic equilibrium; the density distribution in cylindrical coordinates is. 



p{R,z) 



^°(f) 



exp 



R_ 

Rr 



7/2-p-q/2- 



exp 



HiR) 



(1) 



Here H{R) is the scale height at a radius R, which is the sound speed divided by the 
Keplerian angular velocity: H(R) = \/2cs/VL = ^J2kT{R^ 0)R^/GMTh. The surface density 
power-law index can be expressed as 7 = —3/2 + p + q/2, as S(-R) = p{R, 0)H{R)/y/7r. 

The q is the temperature power-law index of dust grains, T{R, 0) = To{Rq/RY. When 
assuming a power-law opacity (k,^ = kq^u / uq)^) and radiative equilibrium with a centra l 
protostar, it is expressed by /3: q = 2/(4 + /3) in the low optical depth limit (ISpitzerlll978l ). 
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For calculating scale heights, we adopted the mid-plane temperature Tm{R,0) = Tq{Rq/R)'^ 
and q = 0.43 (corresponding to (3 = 0.7). Note that the assumed power-law inde x is 
consistent with self-consistent temperature distributions ( iDuUemond fc Dominikll2004al ). In 
contrast to the scale height calculation, we utilized temperature distributions as a function 
of z as well as R, which simulates self-consistent temperatures: 



T{R, z) = WTmiR, 0) + (1 - W)Tsir), 



(2) 



where Tm{R,0) and Ts{r) indicate mid-plane and surface temperature distributions, respec- 
tively. The two temperature distributions have functions of the same power-law index of 
q = 0.43: T^{R,0) = To{Ro/R)i and T,{r) = T,o(r,o/r)«. However, note that the sur- 
face temperature depends on r = y/R"^ + z"^, while the mid-plane temperature is a func- 
tion of R. We adopt T^q = 400 K at r^o = 3 AU corresponding to L = 8.3 L© and 
W = exp[— z^/2(3iJ(i?))^] resulting in temperature distributions close to the self-consistent 
ones empirically, except the very inner region showing a steep temperature gradient. The 
mid-plane temperature distributions (the scale factor To(-Ro)) are most sensitive to the disk 
thickness so they are scaled after checking self -consistent temperature distributions by a 
test run of a specified disk thickness using the iDuUemond fc DominikI ( l2004al ) code. The 
employed values are mentioned in Section 14.11 This estimation of the temperature distribu- 
tion may not be the best approach but it provides a fast means to achieve modeling results 
without a significant loss of accuracy. 

The po of the density expression in Equation ([1]) can be expressed in terms of M^js^ (total 
disk mass), p, q, Rin (disk inner radius), and Re by integrating the de nsity distribution. We 
assume kq = 0.01 cm^g~^ at vq = 230 GHz given by ice mantle grains ([Ossenkopf fc Henning 
19941 ) in a size distribution with a power index of 3.5 ( JMathis et al.lll977l ) and a gas-to-dust 
ratio of 100. The inclination and position angle are also free parameters in our modeling. 



Modeling to constrain these free parameters is carried out with various disk thick- 
ness. For disk thickness we employ a scale height factor, bheight, which is defined as a 
disk scale height with respect to the hydrostatic equilibrium scale height. Therefore, bheight 
smaller than 1 means that the disk appears thinner than the thickness of hydrostatic equi- 
librium in the dust continuum, which implies dust grain settlement. We adopt the cen- 
tral protostellar mass M^yn = 0.55 Mc^ estimated by protostellar evolution tracks on the 



disk gas ( ISargent fc Beckwith 



199ll). 



luminosit y-temperature plot (Beckwith et al.l Il990l ) and by the Keplerian rotation of the 



In summary, the parameters to be constrained for the viscous accretion disk model are 
seven in total: p for a density distribution, (3 for a dust opacity spectral index, Mdisk for a 
disk total mass, R^ for an inner radius. Re for a characteristic radius, 6i for an inclination. 
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and PA for a position angle. In addition, we investigate disk thickness by constraining the 
seven parameters in various hheight cases: hheight = 0.05, 0.1, 0.2, 0.5, 1.0, 1.5, and 2.0. 

Disk images are first calculated by numerically solving the radiative transfer equation 
along the line of sight without optically thin and Rayleigh- Jeans approximations. The result- 
ing model images were "observed" by multiplying by three types of CARMA primary beams 
based on baseline antennas (baselines of 10.1-m and 10.1-m antennas, 10.1-m and 6.4-m 
antennas, and 6.4-m and 6.4-m antennas), Fourier-transforming the image into uv visibility 
space, and sampling with the observational uv coverage. Finally, the sampled visib i lity data 



are c ompared with the observational data in Bayesian inference (JGilks et al.lll996l : iMacKay 



2003|). 



Bayesian inference allows us to obtain the probability distribution of disk properties- 
P{in I D, H), which is a probability distribution of a disk property parameter m with given 
data D and a given disk model H: 

P(,„|^,^).£(iL^^|)M. ,3) 

The P{D I H) called evidence is just a normalization factor of the posterior P{m \ D, H) 
here in our modeling, which employs only a disk model, and we assume uniform pri- 
ors P{m I H) over parameter search ranges. For the likelihood P{D \ m,H), we uti- 
lize Gaussian function, as the noise of interferometric data shows a normal distribution: 
P{D I m,H) = Yl-exp[—{Di — Mi)'^/2af], where Di is an observational visibility. Mi is a 
model visibility, and ai is the uncertainty of the observational visibility. In principle, the 
uncertainty of interferometric visibility data is estimated by system temperatures, band- 
widths, integration tim e, and efficiencies of antenna surface and correlator quantization 



( iThompson et al.ll200ll ). However, these uncertainties do not take into account errors caused 
by atmospheric turbulence, pointing errors, and overall systematic calibration uncertainties. 
We found that the standard deviation of the imaginary components of self-calibrated gain 
calibrator data presents uv distance-dependent uncertainty for the atmospheric turbulence. 
Therefore, we utilize this more-realistic uncertainty, which is typically a few times larger 
than the theoretical estimate. We point out that it is still lower limit of uncertainty, as we 
do not fully capture atmospheric variation for the errors caused in gain solution transfer 
from a calibrator to a target and any other possible error sources, for example, non-identical 
beam patterns of even the same size antennas. 

The parameter search is carried out by a Markov Chain Monte Carlo method (Metropolis- 
Hastings algorithm), and in order to deal with the multi-dimensional parameter space, Gibbs 
sampling is utilized; individual parameters are drawn b ased on the M etropolis-Hastings al- 



gorithm one by one, with the other parameters fixed (JMacKayl l2003l ) . Our modeling has 
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been tested with artificial disk data (lKwonll2009l ). 

In addition to our two wavelength image^j, we take into account the spectral energy 
distribution (SED) of HL Tau over near-infrared to milli meter wavelengths as well, which 



comes from the literature (e.g., iMen'shchikov et al.lll999l ). While geometrical and physical 



properties of a disk are best constrained by high angular resolution images, a best model 
should also explain the SED. The SED is affected by frequency-dependent properties (e.g., 
dust opacity) as well as disk geometries. 



4. RESULTS AND DISCUSSION 



4.1. Disk Thickness : Large Grain Settlement? 

We have obtained marginal probability distributions of the free parameters fitting the 
two images in seven cases with hheight fixed at: bhetght = 0.05, 0.1, 0.2, 0.5, 1.0, 1.5, and 2.0. 
Based on the typical test runs of the seven parameters for each bheight, we adopt Tq = 30, 35, 
41, 55, 70, 82, and 93 K at i?o = 10 AU for the mid-plane temperature distributions. In the 
achieved probability distributions of parameters shown in Figure [21 the case of bheight = 0.05 
is excluded, since the fits always hit our disk mass limit of 0.5 M0. As the bottom-right 
panel of Figure [2] presents log (posterior), the thinner disk models are more preferred by the 
images. The thinner disks have a colder mid-plane temperature so they need more disk 
mass for the observed flux. Note that, however, the disk mass is likely overestimated in thin 
disks, because the gas-to-dust ratio significantly decreases from 100 as dust grains settle on 
the mid plane. As the disk mass goes up, optical depth increases particularly in the inner 
disk region. Therefore, the density gradient gets steeper and the opacity spectral index (3 
becomes larger. The other parameters also change along the disk thickness {bheight) but not 
significantly. 



For the best fit models o: 



each bheight case, we calculat ed model SEDs, using a Monte- 



Carlo radiative transfer code (JDuUemond fc Dominikll2004al ). We assume a power- law opac- 
ity of the constrained index /3 for dust properties. Figure [3] shows the cases of bheight = 0.1, 
0.2, 0.5, 1.0, 1.5, and 2.0 from the bottom, distinguished particularly in the mid-infrared 
regime. The bheight = 1-5 case, which is closest to and does not overesti mate the data points, 



i s pres ented by a solid line. The open rectangular points are data from IMen'shchikov et al. 
(119991 ) and the solid stars mark our values at A = 1.3 mm and A = 2.7 mm. Although the 



^Our original data are visibilities of interferometry and we carried out model fitting in the visibility space. 
However, the word "image" is used to indicate our data in contrast to the spectral energy distribution data. 



thinner disk models are preferred by our two wavelength images of CARMA, they cannot 
explain the mid- infrared fluxes. In contrast, the bheight = 1-5 case reasonably well recovers 
the SED over mid-infrared to millimeter wavelengths and fits the two CARMA images. The 
bump of the model SEDs in near-infrared regime is the cent ral protostar. Since HL Tau is 
located at a boundary of a wall structure (IWelch et al.l 120001 ) . one would expect large-scale 
extinction of those wavelengths. 

The bheight greater than unity can be interpreted as existence of residual envelope com- 
ponents. Indeed, some previ ous studies have repor ted envelope components of HL Tau based 



on SED fitting results (e.g., iRobitaille et al.l 120071 ). However, note that the disk scale height 



larger than the hydrostatic equilibrium value could be an overestimate. Since SED data have 
been taken with poor angular resolution and HL Tau is located on a wall structure, it may 
be likely that the wall structure also contributes to the mid-infrared fluxes. On the other 
hand, it might also be caused by an overestimate of the central protostellar mass of HL Tau 
or an underestimate of the disk mid-plane temperature. 

It is noteworthy that thin disk models are preferred by the image fitting, while thick 
disk models better match SED. We do not attempt to quantitatively fit the two images and 
the SED simultaneously in Bayesian inference, as it is not clear how the significance of data 
should be distributed between the images and the SED, which consist of millions of visibility 
data points and tens of flux densities, respectively. More importantly, protoplanetary disk 
emission is unlikely to be explained by a simple model with a scale height. In other words, 
disk thickness {bheight in our modeling) can depend on dust grain sizes, as larger grains likely 
settle to the mid plane faster and more compactly. This results in smaller scale heights 
for longer wavelength observations; short wavelength observations do not probe the mid 
plane due to relatively large optical depth, so a thick disk would be preferred. In fact, 
numerical simulations show that large grains rapidly settle i nto the disk mid plane, resulting 



in sig nificant differences of scale heights along grain sizes ( jBalsara et al.l 120091 : iTilley et al. 



20101 ) . This is not surprising since for a l ong time dust settlement has been considered as the 



initial stage to form planetesimals (e.g., iGoldreich fc Wardlll973[ ). We also found the best 



fitting model of only our millimeter images with setting the bheight as a free parameter and 
with a broader disk mass range. When the mid-plane temperature is assumed as Tq = 30 K 
at 10 AU, which is the same to the bheight = 0.05 case, the marginal probability of bheight is 
at maximum around 0.02 and the other parameters follow the trend shown in Figure |2l The 
results imply that the dust grains most sensitive to millimeter wavelength observations have 
settled in a scale height deceased by a factor of 50. 



It has been noted that dust settlement decreases mid- and far-infrared emission (JDuUemond fc Dominik 



2004bl ). as shown in the thin disk cases of Figure [31 However, HL Tau has bright mid-infrared 
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emission; a disk that still has much gas and "strong" accretion like HL Tau, which has rela- 
tively strong turbulence, can arguabl y keep significant ra id-infrared emission by not-settling 



fine dust grains mixed with gas (e.g., iBalsara et al.ll2009[ ). On the other hand, we are aware 



that the disk thickness can also depend on disk models (e.g., accretion disk and power-law 
disk), so it may be worthy of testing different disk models further (Kwon et al. in prepara- 
tion). 



4.2. Disk Properties Constrained by Modeling 

Considering the fitting results of the two wavelength images and SED, the possibility of 
large grain settlement, and the bias of the extended wall structure, it is not trivial to select 
the best fitting model. In addition, if the disk is stratified as discussed, it might not be the 
best approach either. If we simply took account of the fitting results of the two images and 
SED, however, the bheight = 1-5 would be a representative of the best fitting models, with a 
caveat that the disk could be stratified with large grains settled down on the mid plane. In 
addition, as dust is more or less 1% of the total disk mass, the bheight = 1-5 case may better 
present the overall disk structure of gas and dust. We discuss each of the model parameters 
focusing on the bheight = 1-5 case in the following paragraphs. On the other hand, the 
stratified structure may be investigated using submillimeter observations with high angular 
resolution and image fidelity, e.g., using the Atacama Large Millimeter Array (ALMA). Our 
A = 2.7 mm image does not have a high angular resolution comparable to A = 1.3 mm. If 
HL Tau is apparently stratified, a thicker disk than the one constrained by our millimeter 
data would be preferred by submillimeter data. 

Figure H] presents uv visibilities in annulus average with the dust opacity spectral index 
P simply calculated assuming optically thin and Ray leigh- Jeans approximation: /3 = a — 2, 
where a is a fiux density spectral index. The overlaid solid lines are of the best fit model with 
bheight = 1-5. Figure [5] shows the image, model, and the residuals. The residuals to the model 
at A = 1.3 mm are a mixture of plus-minus 3cr peaks; the residuals to the A = 2.7 mm fit 
show a 4a peak to the northeast and a 7a peak to the southwest. The A = 2.7 mm residuals, 
which are formally significant, might indicate dust opacity variation in the disk, but it is not 
clear due to our poor angular resolution in the A = 2.7 mm map. On the other hand, the 
residuals could be free-free emissio n from the ionized jet oriented at this position angle and 



seen better at longer wavelengths (JRodmann et al.l 120061 ). 



The probability distributions for the free parameters constrained by our modeling with 
bheight = 1-5 are presented in Figure [6] and in Table [TJ The parameter ranges for good fits in 
the figure and the uncertainties in the table are for our specific parameterized accretion disk 
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models; other families of models may also fit the data. In addition, the uncertainties of the 
figure do not include systematic uncertainties such as absolute flux calibration uncertainty. 
The table has a column and notes for possible systematic uncertainties. 

The disk central density power-law index p is well constrained in the model at p = 1.06, 
under the temperature power-law index, q = 0.43. Based on the fitted density power-law 
index, a surface density distribution power-law index 7 of -0.22 is derived. This negative 
value does not imply a continuous increase of surface density with radius, due to the tapering 
of the exponential function. For the best fit model in Table [H the surface density increases 
by 50% from Rin to 28 AU, then decreases beyond 28 AU. While the density distribution 
is not straightforward to be compared with previous studies using a power-law disk model 
due t o the different functio nal forms, a study using the accretion disk model toward HL 



Tau ( iKitamura et al.l 120021 ) reported a steeper surface density gradient (7 = 0.78). This 



discrepancy may be caused by the poor angular resolution of their observations, 1'.'2 x I'.'l. 

The dust opacity spectral index /3 is best fitted by a value of 0.73. This value is primarily 
determined by the ratio of the flux measure at the two wavelengths since our model fits do 
not have significant optical depths over most of the disk. As shown in Figure HI however, 
a thinner disk has significant effects of optical depths resulting in a larger /3. The total 
uncertainty in /3 includes the statistical uncertainty, which is captured in the figure, and the 
relative uncertainties in the flux calibration between the two wavelengths. For the estimates 
flux uncertainties of 8-10% for the two wavelengths, the total uncertainty in /3 is ibO.25. 
The small (3 implies that dust grains have grown furthe r comparing protostellar systems at 
earlier stages, which have /3 ~ 1 (e.g. jKwon et al.ll2009l ). 



On the other hand, iRodmann et al.l (120061 ) derived (3=1.3 in the wavelength range of 1 



to 7 mm. Extension of our /3=0.73 law to 7 mm clearly overestimate s the measured f l ux: 9. 4 



mJy from our model versus 3.9 mJy of disk emission estimated by iRodmann et al.l (120061 ). 
To produce the measured 7 mm flux our model needs /3 = 1.28 between A = 1.3 mm and 
7 mm. Note that the model with this /3 produces the same A = 1.3 mm image but it does 
not fit our A = 2.7 mm flux and image. Figure [7] shows the ratio of the mass absorption 
coefficient {k,^) at a wavelength divided by the mass absorption coefficient at A = 1.3 mm. A 
straight line in this log-log plot represents a single power law index over the full wavelength 
range. The turnover in th e ratio value could indicate a largest grain size of several millimeter 



in the disk (JDraind 120061 ). 



The disk mass is very well determined within the context of the model: 0.135 Mq. How- 
ever, the full uncertainty must include uncertainty in the absolute flux calibration, the mass 
absorption coefficient, the gas-to-dust ratio, and the dust temperature. We have estimated 
the flux calibration uncertainty at 10% which translates directly to a mass uncertainty. The 



- 11 - 



mass opacity coefficient could be uncertain at the factor of two level. The total luminosity 
of the HL Tau system is also uncertain: estimates range from 1 to 11 Lq. The resulting 
uncertainty in the dust temperature at a fiducial radius is about 67%. Note that the massive 
disks constrained in thinner disk models (Figure [2]) would have larger uncertainty in the 
gas-to-dust ratio, as the models imply grain settlement. Although the mass is uncertain, 
we argue that this particular model of hheight = 1-5 would have the least overestimated disk 
mass and the least uncertainty in our modeling at about a factor of two. 

The inner radius is constrained by our data to be significantly smaller than our A = 1.3 mm 
linear resolution of about 20 AU. The formal model fit gives an inner radius of 2.4 AU. Fits to 
the mid- infrared emission typically require disk material to an inner radius of several stellar 
radii to 0.2 AU ( iMen'shchikov et al.lll999t iFurlan et al.l 120081 ). It is possible that an inner 
radius in millimeter emission is tracing a surface density drop, but the material remains 
optically thick at mid-infrared wavelengths throughout the inner disk. 

The characteristic radius, which divides the regions dominant by the power-law and 
the exponential density distributions, is 79 AU. It is this exponential cutoff in Equation 
([T]) that sets the outer scale for the disk in the presence of the shallow power-law index. 
The surface density falls by a factor of ten from its peak at a distance of 125 AU, and a 
factor of one-hundred at 165 AU. In the context of the accretion disk model, the transitional 
radius Rt-, the radius where the bulk flow direction changes from inward to outward for 
angular momeiitum c onservation, is then 40.3 AU: Rt = Rc{S/2Y, where 6 = 1/(2 — 7) 



flAndrews et al.l 120091 ). A previous study reported that Rt gets larger with age fllsella et al. 



20091 ) but another circumstellar disk survey did not find such a trend ( JAndrews et al.ll2009l ). 



The disk inclination and position angle are very well determined by the disk model: 
inclination of about 40° (face-on disk 6i = 0°) and position angle of about 136°. In fact, 
these are not varying much even over the wide range of hheight (Figure |2]). The inclination is 



particularly consistent with previous studi es using submi l 



observations of a hi gh angular resolu tion. iMundy et al.l (Il996() suggested a rela t ively wide 



imeter and millimeter wavelength 



range of 20°-55° and lLav et al.l fll997f ) constrained as 42° ±5°. IMen'shchikov et al.l f ll999l ) also 
reported an inclination of 47° using a detailed SED modeling. However, fits at ir ifrared wave- 
lengt hs have suggested inclinations from i = 15° from mid-infrared spectr al fits (iFurlan et al. 



20081 ). to 66-71° from fitting to near- infrared images and polarization ( iLucas et al.l 120041 ). 
It is possible that the near- and mid-infrared determinations of inclination are affected by 
large scale inhomogeneities such as the wall structure. 
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4.3. Constraints on the Presence of the Proposed Protoplanet and Disk 

Stabihty 



Figure IHpresents the zoomed-in residual map overlaid on top of the image at A = 1.3 mm. 
Particularly the top-left panel displays a blow-up of the A = 1.3 mm image of the disk in 
the region of the protoplanet candidate (JGreaves et al.l 120081 ) . The right panel shows the 
residuals for the same region when t he disk model is rern oved. The cross-hashed circles of 
both panels rnark th e position where iGreaves et al.l ( l2008l ) claimed a protoplanet candidate. 
Greaves et al.l ( l2008l ) detected a clump of 78 ± 17 fiJj at A = 1.3 cm. The clump would be 
25 ±5 mJy at A = 1.3 mm, assuming a z/^'^ emission spectrum. If the dust grains responsible 
for the emission were very large so that the emission spectrum rather followed i/^°, i.e., the 
hardest to detect for our observations, then the flux at A = 1.3 mm would be 7.8 ± 1.7 mJy. 
Since our observational sensitivity at A = 1.3 mm is 0.8 mJy beam~^, the candidate could be 
detected in our observations in spite of our relatively worse angular resolution: 0.08" versus 
0.13". 

As shown in the panels, the candidate position is within the dust thermal emission region 
of the disk at A = 1.3 mm; however, subtraction of the axi-symmetric disk model removes 
essentially all flux. Before model subtraction the flux at the position is ~16 mJy beam~^; 
after subtraction the position has a negative flux at a level of 2a. Therefore, our A = 1.3 mm 
image provides no support for the candidate protoplanetary condensation. The flux level 
before subtraction is roughly consistent with the value bootstrapped from the A = 1.3 cm 
flux, but it is part of the large disk structure, which was not detected by their observations, 
rather than a compact signal. 

Although our observations do not have the protoplanet candidate signal, the HL Tau 
disk could be gravitationally unstable. Our disk model provides a surface density profile that 
can be used to calculate local gravitational stability in the disk. The Toomre Q parameter 
( lToomrelll964[ l with values smaller than 1.5 is a good indicator for gravitationally unstable 
regions. The Q parameter versus radius. Figure IHl was calculated with the disk temperature 
and density distributions of our modeling. It shows a minimum less than or close to 1.5 
in the range of 50 to 100 AU, which means that substructures can dev elop at these radii 
by gra vitational instability. This result is in agreement with the work by iNero &: Bjorkman 
( 120091 ). which studied the possibility of fragmentation in the HL Tau disk based on perturba- 
tion cooling time and gravitational instability. Considering the surface density, protostellar 
mass, and the radius of our modeling in the iNero fc Bjorkmanl (120091 ) results, roughly 2 
Mjupiter fragments are favored between 50 and 100 AU. As discussed in Section 14. 2[ the 
disk mass that is the most important property for the Q parameter calculation has a large 
uncertainty. However, since the model used here is not even the one with a small hheight 
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favoring a massive disk, the HL Tau disk is likely unstable in the outer region. 

Although we have not detected a clear substructure, the HL Tau disk is likely unstable 
and can be fragmented. The positive residual at ~100 AU in Figure [S] niight be a hint of a 
dust lane. I ndeed, extrasolar planetary systems w ith planets at outer disk regions have been 



found (e.g., iKalas et al.ll2008l : iMarois et al.ll2008l ). These planets of outer disk regions are 



thought to have formed by gravitational instability (so-called the disk instability scenario), 
although it is little known whether and how the protoplanets formed by gravitational insta- 
bility at the early phase of protoplanetary disks s till having gas can survive the following 



active migration stage (e.g., iD'Angelo et al.ll2010l ). In contrast, the other major model of 



giant planet formation, so-called the core accretion scenario, cannot currently explain the 
giant planets at the outer disk regions due to too-large time scales of gas accretion (e.g.. 



D'Angelo et al.ll2010l ). Regardless of giant planets being actually formed by disk instability 
or core accretion, they should be formed in protoplanetary disks that still have gas. In the 
near future, ALMA will provide the unprecedented image to reveal the more detailed disk 
structure of HL Tau. In addition, we would be able to distinguish the two scenarios of giant 
planet formation by gas kinematics and substructure and scale height variation in multiple 
wavelength observations. 



5. CONCLUSION 

We have obtained dust continuum data of HL Tau at A = 1.3 mm and 2.7 mm with 
the exceptionally high image fidelity and angular resolution up to 0.13" using CARMA. 
This resolution reveals 18 AU scales at the distance of HL Tau, 140 pc. Through visibility 
modeling in Bayesian inference adopting the viscous accretion disk model, we constrained 
the physical properties of HL Tau. As we fit our two wavelength data simultaneously, we 
also constrained the dust opacity spectral index. In addition to the two wavelength images, 
we also qualitatively fit the SED. 

Modeling in Bayesian inference toward the two wavelength images prefers a thin disk 
(Figure [2]). However, thin disk models can not explain mid- infrared emission of SED (Figure 
|3]). A thick disk is needed to produce the mid-infrared emission. The discrepancy between the 
fitting results of millimeter images and the SED may indicate the settlement of particularly 
large grains into the mid plane (a stratified settlement of dust grains). However, based on 
the reasonable fitting to both the millimeter images and SED and considering that dust is 
not the majority of disk mass, we selected and further discussed the case of hheight = 1-5 
(Figure |6] and Tabled]), which possibly supports envelope residuals beyond the hydrostatic 
equilibrium thickness. The scale height ratio larger than unity might also be caused by an 
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overestimate of the HL Tau protostellar mass or an underestimate of the disk mid-plane 
temperature. The bheight = 1-5 model reveals that the disk mass is about 0.135 Mq and the 
dust opacity spectral index presents that grains have grown further, compared with objects 
at earlier evolutionary stages. The disk inclination and position angle are well determined 
as 40° (face-on disk of 6i = 0°) and 136°, respectively. 

Although a protoplanet candidate has been claimed in HL Tau, our observations, which 
are more sensitive to dust thermal emission, do not detect emission structure supporting 
the candidate. However, HL Tau appears to be gravitationally unstable and could favor 
fragmentation on the Jupiter mass scale between about 50 and 100 AU. Further observations 
with a higher angular resolution and a higher sensitivity using ALMA will be able to show 
substructures possibly formed by gravitational instability as well as to verify the stratified 
settlement of dust grains in HL Tau. 
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Fig. 1. — HL Tau in the A = 1.3 mm continuum. The image is the combined of CARMA A, 
B, and C configurations, and the synthesized beam is 0'.'17 x 0'.'13 {PA = 86°) corresponding 
to 18 AU. The contour levels are 2.5, 4.0, 6.3, 10, 16, 25, and 40 times a = ±0.8 mJy beam~5 
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Fig. 2. — Free parameter variation along changing bheight- The thinner models are preferred 
by our two wavelength images, based on the posterior values. However, thin disk models 
cannot explain the mid-infrared fluxes (Figure [3]). 
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Fig. 3. — HL Tau SED overlaid with models of various hheight values. The solid line is the 
case of hheight = 1-5 and the dashed lines are cases of hheight = 0.1, 0.2, 0.5, 1.0, and 2.0 from 
the bottom. The two stars indicate our data points 700 and 120 mJy at A = 1.3 and 2.7 
mm, respectively. 
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— HL Tau visibility averaged in annuli. The rectangular and triangle points are data 
1.3 mm and 2.7 mm, respectively, and the solid line is the model best fitting to both 
and SED. 
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Fig. 5. — HL Tau continuum, model, and residual maps at A = 1.3 mm and 2.7 mm. The 
A = 1.3 mm image is the same as Figure [1] and the synthesized beam of the A = 2.7 mm 
image is 0l!98 x 0'.'70 {PA = 80°). The contour levels are 2.5, 4.0, 6.3, 10, 16, 25, and 40 
times a = ±0.8 and ±1.1 mJy beam^^ at A = 1.3 mm and 2.7 mm, respectively. 
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Fig. 6. — Posterior distributions of parameters in HL Tau. The solid lines of the plots 
indicate normal distributions of the posterior-weighted means and standard deviations of 
the parameter distributions. The posterior-weighted means and standard deviations are 
listed in Table [1] with the parameter search ranges. 
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Fig. 7. — Mass absorption coefficients [ku) to produce the measured fluxes at lower frequen- 
cies with our disk model. The slopes present dust opacity spectral indexes (/3) determined 
between a frequenc y and v = 230 GHz. Th e circles mark ou r resu lts, the triangles are of 
Wilner et al.l ( 1l996l ). and the squares are of iRodmann et al.l (120061 ). The open and ffiled 
symbols present the mass absorpti on coefficients when free-free emission components are in- 
cluded and excluded, respectively. IWilner et al.l (119961) estimated 3 0% of their flux 10 ± 1.4 



Rodmann et al.l (l2006l ) argued 20% of their flux 



mJy at A = 7 mm for free-free emission and 

4.92 mJy at A = 7 mm and 50% of 1.63 mJy at A = 1.3 cm for free-free emission. 
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Fig. 8. — HL Tau in the A = 1.3 mm continuum overlaid with residual contours in the 
bottom. The disk image is the same as Figure [1] and the levels of residual contours are 2 
and 3 times a = ±0.8 mJy beam~^. The top row shows the northwest quadrant region 
zoomed-in. The contour levels of the left are 2.5, 3.0, 4.0, 5.0, 6.3, 8.0, 10.0, 11.5, 13.0, 14.5, 
16.0, 18.0, 20.0, 22.0, 25.0, 28.0, 32.0, 36.0, and 40.0 times a = ±0.8 mJy beam^^ The red 
hashed circles indicate where a protoplanet has been claimed. As shown, no excess signal 
has been detected. 
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HL Tau: accretion disk model 
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Fig. 9. — Toomre Q parameter along radius in HL Tau. The dotted lines indicate the region 
when assuming factor of two uncertainty in the disk mass mainly induced by kq uncertainty 
( lOssenkopf fc Henninglll994l ). Note that the Q value is smaller than or very close to 1.5 
between 50 and 100 AU, which suggests that the region may be gravitationally unstable. 
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Table 1: Disk fitting results of free parameters. 



Parameters 


Search ranges 


Fitting results 


Statistical 


Systematic 










Uncertainties 


P 




0.00-3.35 


1.064 


0.0085 


0.07^ 


P 




-0.35-1.7 


0.729 


0.0069 


0.25*^ 


Mdisk 


[Mol 


0.005-0.2 


0.1349 


0.00034 


0.013= 


Rin 


[AU] 


0.1-20 


2.4 


0.45 


2.1^^ 


tic. 


[AU] 


40-460 


78.9 


0.20 


9^ 


Oi 


o 


0-83 


40.3 


0.23 




PA 


o' 


0-180 


135.8 


0.35 




bheight 




0.05-2.0 


1.5 


fixed 


7 






-0.22 


derived from 


Rt 


[Au; 




40.3 


fitted parameters 















"Roughly estimated assuming that A-configuration data have a 10% flux cahbration error with respect to the 
other configuration data. Note that relative errors in flux calibration over different configurations can cause 
a gradient change along radius in brightness. The p is dependent on the assumed functional description in 
Equation ([T]) and the assumed value of q. 

^Estimated from the assumed 10% and 8% uncertainties in the A = 1.3 mm and A = 2.7 mm fluxes, respec- 
tively. 

'^Estimated using the assumed flux calibration uncertainty in the A = 1.3 mm flux. Note that the disk mass 
is sensitive to the temperature scale (Tq), the dust mass absorption coefficient, and the gas-to-dust ratio. 
Therefore, the total uncertainty in the value can be up to about a factor of two. 

''Estimated from mapped model images with a varying inner edge of the disk to give a detectable emission 
decrease in the central position. 

'^Estimated from half of the angular resolution. 



